library(boot)
library(readr)
library(pscl)


moore <- read_csv("moore_toreplicate.csv")

civil_war_form <- refugees ~ civwar
moore_form <- refugees ~ deathmag + violdissent+ civwar+ waronter +dem_aut+ trans+ bothgnp + ptssd 
our_form <- refugees ~ deathmag +violdissent+ civwar+ waronter +dem_aut+ trans+ bothgnp + ptssd + gdppc_weighted 
our_polity_form <- refugees ~ deathmag +violdissent+ civwar+ waronter +dem_aut+ trans+ bothgnp + ptssd +  polity_weighted 

countries <- unique(moore$ccode)
actual <- moore_output <- our_output <- civil_war_output <- our_polity_output <- NULL
seed <- 123
for(i in 1:length(countries)){
  if(i%%10 == 0) {print(i)}
  learn <- test <- NULL
  samp <- countries[i]
  learn <- moore[!moore$ccode %in% samp, ]
  test <- moore[moore$ccode %in% samp, ]          
  lm.civil.war <- lm.moore <- lm.our <- lm.our.polity <-  NULL
  
  try(lm.civil.war <- zeroinfl(civil_war_form, data=learn, dist='negbin'))
  try(lm.moore <- zeroinfl(moore_form, data=learn, dist='negbin'))
  try(lm.our <- zeroinfl(our_form, data=learn, dist='negbin'))
  try(lm.our.polity <- zeroinfl(our_polity_form, data=learn, dist='negbin'))
  
  try(test$pred.civil.war <- predict(lm.civil.war, type='response', newdata = test))
  try(test$pred.moore <- predict(lm.moore, type='response', newdata = test))
  try(test$pred.our <- predict(lm.our, type='response', newdata = test))
  try(test$pred.our.polity <- predict(lm.our.polity, type='response', newdata = test))
  
  actual <- c(actual, test$refugees)
  civil_war_output <- c(civil_war_output, test$pred.civil.war)
  moore_output <- c(moore_output, test$pred.moore)
  our_output <- c(our_output, test$pred.our)
  our_polity_output <- c(our_polity_output, test$pred.our.polity)
}

results.civil.war <- abs(civil_war_output - actual)
results.moore <- abs(moore_output - actual)
results.our.polity <- abs(our_polity_output - actual)
results.our <- abs(our_output - actual)

par(mar=c(13,3,2,1), pty='s')
b1<-barplot(c(median(results.civil.war, na.rm=T),
              median(results.moore, na.rm = T),
              median(results.our.polity, na.rm = T),
              median(results.our, na.rm = T)), 
            main = '', las=1,   xpd=F, col='white') 



axis(1, at=c(0.2, 7.2), labels=FALSE, lwd.ticks=0)
axis(1, labels = c("Base","Moore & S","Our Polity","Our GDPPC"),
     las = 2, cex.axis = 1.7, at = b1)
mtext('Median Absolute Error', 2, line=4, cex=1.7)
boot.median <- function(data, indices){
  d <- data[indices] # allows boot to select sample
  return(median(d))
}
# bootstrapping the CI of the median     
b1.index <- 0
ups <- lows <- NULL
for(this.data in list(results.civil.war[!is.na(results.civil.war)],
                      results.moore[!is.na(results.moore)],
                      results.our.polity[!is.na(results.our.polity)],
                      results.our[!is.na(results.our)]
)){          
  b1.index <- b1.index+1
  print(b1.index)
  results <- boot(data=this.data, statistic=boot.median, R=5000)
  # get 95% confidence interval 
  ci.95.low <- unlist(boot.ci(results, type="bca"))$bca4
  ci.95.high <- unlist(boot.ci(results, type="bca"))$bca5
  ups <- c(ups, ci.95.high)
  lows <- c(lows, ci.95.low)
  # add a line for the 95% CI (and center it)
  CI.size <- ci.95.high - ci.95.low
  segments(b1[b1.index], median(this.data[!is.na(this.data)]) + CI.size/2, b1[b1.index], median(this.data[!is.na(this.data)]) - CI.size/2)
}  
